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Abstract 

Bezier curves are a widespread tool for the design of curves in Euclidian space. 
This paper generalizes the notion of Bezier curves to the infinite-dimensional space 
of images. To this end the space of images is equipped with a Riemannian metric 
which measures the cost of image transport and intensity variation in the sense 
of the metamorphosis model IMY01I . Bezier curves are then computed via the 
Riemannian version of de Casteljau’s algorithm, which is based on a hierarchical 
scheme of convex combination along geodesic curves. Geodesics are approxi¬ 
mated using a variational discretization of the Riemannian path energy. This leads 
to a generalized de Casteljau method to compute suitable discrete Bezier curves 
in image space. Selected test cases demonstrate qualitative properties of the ap¬ 
proach. Furthermore, a Bezier approach for the modulation of face interpolation 
and shape animation via image sketches is presented. 


1 Introduction 

Bezier curves are a classical tool in computer aided geometric design. We generalize 
this tool to an infinite-dimensional manifold of images. To this end we consider images 
as control points and compute a curve in the space of images by a Riemannian version 
of de Casteljau’s algorithm. Thus, our approach relies on the definition of a Riemannian 
structure and the resulting notion of geodesic curves on the space of images. 

The study of the space of images from a Riemannian manifold perspective allows 
to transfer various tools from classical geometry and computer-aided design to this 
infinite-dimensional space. E.g., geodesics as minimizers of the path length generalize 
the notion of straight lines from Euclidian space, the geodesic distance allows to mea¬ 
sure dissimilarity of images in a rigorous way, the geometric exponential map allows 
to generate natural extrapolation paths starting from an infinitesimal variation of an 
image. During the past decade, this geometric approach triggered the development of 
new methods in computer vision and imaging, ranging from shape statistics 1FLPJQ41 
to computational anatomy I BMTY02 1. 

Different image manifolds have been investigated. The concept of optimal transport 
was used to study the space of images, where image intensity functions are considered 
as probability measures I B BOO , ZYHT07 1. Following the classical paradigm by Arnold 
I Arn66l 1AK98IL the temporal change of images can be studied based on the flow of 
diffeomorphism concept, where a family of diffeomorphisms (^(t))te[o,i] : D 
R d on D C describes a flow transporting image intensities along particle paths. 
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Figure 1: Given 4 portrait images as control points (A-D) a discrete Bezier curves in 
the metamorphosis manifold (marked in blue) are compared with the geodesic interpo¬ 
lation (marked in green). 


The metamorphosis approach was first proposed by Miller and Younes 1MY011 as a 
generalization of the flow of diffeomorphism model and comprehensively analyzed by 
Trouve and Younes HTY05bl . Besides pure transport it allows in addition for image 
intensity variations along motion paths. 

On any Riemannian image manifold geodesics are not only minimizers of the 
length functional but also of the path energy, which is the time integral of the squared 
path velocity. In addition, minimizers of the path energy have constant speed param¬ 
eterization. Hence, these geodesics are not only the obvious generalization of straight 
lines in Euclidian space, but also allow a simple procedure to compute convex combi¬ 
nations with convex coefficients t and 1 — t: one simply evaluates the path in image 
space at time t£, where £ is the total path length. Convex combination is the key ingre¬ 
dient of de Casteljau’s algorithm. This is our starting point, and in this paper we will 
focus on the metamorphosis model and use the above insight to study Bezier curves in 
this space. 

Already in 1995, Park and Ravani realized how de Casteljau’s classical algorithm 
can be applied on Riemannian manifolds (just replacing straight lines by geodesics) to 
define smooth Bezier curves, and they applied this concept to generate smooth trajec¬ 
tories for moving objects on a manifold of kinematically admissible motions 1PR951 . 
Such Bezier curve segments on manifolds can also be patched together to yield a spline 
interpolation of points on the manifold. To this end, Popiel and Noakes derived formu¬ 
lae for the endpoint velocities and (covariant) accelerations of Bezier curves in terms of 
the Bezier control points I PN07b L which can be employed to ensure smoothness of up 
to two derivatives. In IMCV08 L a concept of geodesics from discrete differential ge¬ 
ometry is used to apply de Casteljau’s algorithm for the purpose of curve modeling on 
triangulated surfaces. C 1 smoothness of Bezier splines is here ensured by extrapolat¬ 
ing the endpoint velocity of a Bezier patch via so-called straightest geodesics to obtain 
control points of the next patch. The authors also implement a classical control point 
refinement procedure to approximate Bezier curves via a subdivision scheme. Gousen- 
bourger et al. compute interpolating C 1 splines on the sphere, SO( 3), and a shape 
manifold of closed curves 1GSA141 . where they generalize the idea of cubic Bezier 
splines with minimum acceleration to manifolds. Let us finally mention that there are 
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also alternative generalizations of polynomial curves to manifolds based on variational 
definitions. In particular, generalized cubic curves can be defined as minimizers of the 
integrated squared second covariant derivative I CSC0 1, PN07al . 

In this paper we discuss Bezier curves in the space of images considered as a Rie- 
mannian manifold using the metamorphosis approach. Figure [T] shows the comparison 
of piecewise geodesic interpolation and a discrete cubic Bezier curve for four different 
facial images as control points—from now on called control images. The Bezier curve 
is temporally smooth, comes close to the two intermediate images, and features of the 
control images globally pervade the curve. 

2 Time Discrete Metamorphosis 

In this section, we briefly review the metamorphosis model and introduce a suitable 
variational time discretization. With the notion of continuous and discrete path energy 
at hand we can define a continuous and discrete geodesic interpolation. 

2.1 Continuous Model 

The Riemannian metric in the metamorphosis model combines a measurement of fric¬ 
tion caused by the flow of the image intensities and a measure of the variation of in¬ 
tensity values along motion paths. More explicitly, along a path u : [0,1] -> L 2 (D) in 
the manifold of images on a bounded Lipschitz domain D we consider the Riemannian 
metric 2 

9 u(t){u{t),u(t)) = min j^L[v(t),v(t)] + i da;, (1) 

where ii(t) represents the velocity along the path (and is just the pointwise time deriva¬ 
tive of the image u{t)). The first term describes the viscous dissipation in a multipolar 
fluid model (c/. EMD)- Here, we assume that velocities field vanish on dD and 
consider L[v , v\ := ^(ti£[v]) 2 + fiti{e[v] 2 ) + 'y\D m v\ 2 as the classical viscous dissi¬ 
pation model for a Newtonian flow plus a simple multipolar dissipation model, where 
e[v] = §(Vu + V?; T ), m > 1 + § (for d the space dimension) and A, /i, 7 > 0 . 
The second term with weight ^ > 0 measures the temporal variation of the intensity 
in terms of the material derivative ^ u = ii -\-\7u • v. Obviously, the same temporal 
change u(t) in the image intensity can be produced by different motion fields v(t) and 
associated material derivatives which makes the minimization with respect to v in 
the definition of the metric necessary. The path energy is given by 

£[(^M)te[o,i]]= [ 9u{t)(u(t),u(t))dt . ( 2 ) 

Jo 

Dupuis et al. I DGM98 1 showed already for the pure transport model (5 = oo) that 
paths of finite energy are indeed one-parameter families of diffeomorphisms. A rigor¬ 
ous analytical treatment of the general model can be found in [ TY()5a l including the 
existence of minimizing paths connecting two images in L 2 (D). It relies on an appro¬ 
priate notion of tangent vectors, which are equivalence classes of pairs (u, z), where 
v is a motion field with square integrable Jacobian, and 2 is a weak representation of 
the material derivative m u ■ Let (*u(tj |q denote a continuous geodesic curve in the 
space of images given as the minimizer of the path energy for given images iz( 0 ) = ua 
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and u( 1) = ub , then we can define a convex combination of the images ua and ub 
with weight A via an interpolation of the geodesic at time t = X,i.e. 

l(u A ,u B ,X) =* «(A). (3) 

Here, we in particular make use of the fact that the “speed” \Jo u (i) (Mt)- W) ) is con- 
stant for a minimizer of the path energy. 

2.2 Time Discrete Model 

To robustly and efficiently approximate geodesic curves in the metamorphosis model 
we define a suitable time discrete version of the continuous path energy To this 
end, we consider (K + 1)-tuples (uo, ..., uk ) of images and define on them an energy 

K 

E[(uo, ■ ■ ■ , u K )] as K w[u k - 1 , Uk] , (4) 

A—1 

where W^-i, is a classical matching functional measuring the cost to match the 
image Uk-i with the image Uk- In detail, 

W[u,u] = min [ W(D<j)) + ')\D m (j)\ 2 + \\u o cj) — u\ 2 dx (5) 

<t< Jd ° ' 

for an isotropic and rigid body motion invariant energy density W. A suitable choice 
in the case d = 2 is given in I BER15 I. In our numerical computations we will use 
the simplified energy W{D(j>) = \D<p\ 2 and we will replace 7|Z) m 0| 2 by y|A0| 2 
corresponding to the quadratic form L[v(t),v(t)] = Dv : Dv + ^Av • Av. Mini¬ 
mization is performed over a set of admissible deformations with <fi(x) = x on dD. 
A discrete geodesic path is defined as a minimizer of the discrete path energy E for 
uq = ua and uk = ub- In this discrete path energy, two opposing effects can be 
observed. The last term in ([5]) penalizes intensity variations along the discrete motion 
path (x, (x), (02 ° 0i) , (0 k o...o0 1 )(x)), whereas the first two terms pe¬ 

nalize deviations of the (discrete) flow along these discrete motion paths from rigid 
body motions. In fact, K(uk o 0& — Uk-i) plays the role of a time discrete material 
derivative along the above discrete motion path. This ansatz ensures that the discrete 
energy E T-converges for K —oc to the continuous energy £. For the proof we 
refer to I BER1 51. In particular, discrete geodesics converge to continuous geodesics 
for K oo. An introduction to variational time discretization on shape manifolds can 
be found in I RW14 1. Obviously, the energy scales quadratically in the deformation 0, 
which itself is expected to scale linearly in the time step r = . This motivates the co¬ 
efficient K in front of the discrete path energy. Now, assuming uniqueness of discrete 
geodesics a discrete geodesic interpolation in analogy to the continuous interpolation 
is defined as 

I K (u A ,u B ,k) = u k (6) 

for the discrete geodesic path (uo, ..., uk)- 

3 De Castelj au’s Algorithm on the Manifold of Images 

A Bezier curve is a polynomial curve of degree n > 1 in Euclidean space R d that 
is defined by n + 1 control points xq, ..., x n E M d . The curve emanates from the 
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first and ends in the last control point, but does not interpolate the intermediate control 
points. However, it roughly follows those points and always remains within their con¬ 
vex hull. From another viewpoint, the choice of the n + 1 control points is equivalent 
to specifying the first derivatives at the curve end points. The classical de Castel- 
jau algorithm for evaluating points on a Bezier curve can be used to generalize those 
curves from M n to Riemannian manifolds 1PR951 . 

3.1 Continuous Bezier curves 

De Casteljau’s algorithm constructs each point along a Bezier curve recursively via 
an iteration of weighted interpolations between point pairs, starting from weighted 
interpolations between all neighboring control points. Since weighted interpolation 
between points can also be performed on Riemannian manifolds, the algorithm can 
directly be applied to control points on a manifold, yielding a generalisation of Bezier 
curves. Note that on a Riemannian manifold, weighted interpolation between two 
points yields a point along the connecting geodesic as opposed to a point along the 
connecting straight line in Euclidean space. 

In our case, the control points are images Uq,...,u^ : D R, and the correspond¬ 
ing Bezier curve t B(uq, ..., t) at any time t G [0,1] is recursively defined via 
de Casteljau’s algorithm as 


B(ui, .. .. 


li t),B(u i + li ...,u j ,t),t). 


Geodesics are not unique in general, hence the interpolations and the resulting Bezier 
curve for a given set of input images are in general non-unique as well. We obtain the 
following existence and continuity result. 

Theorem 3.1 (Existence and stability of Bezier curves). Let D be a Lipschitz domain 
and n > 1. For any n + 1 input images Uq, ..., G L 2 (D) there exists a Bezier 
curve B(uq, ... •) : [0,1] L 2 (D). Furthermore, if(u ^ k ,..., u^ k )k=i,2,... is a 

sequence of control points converging against (v ,®,..., u^) in (L 2 (D)) n+1 , then there 
exist Bezier curves B(u [] ..., fc , •) of which a subsequence converges pointwise 
against a Bezier curve B(uq, ..., iz® , •). 

Remark 3.1.1. In the case that all involved geodesics and thus all Bezier curves are 
unique, we even have convergence of the full sequence. 

Proof The result automatically follows from the recursive definition of Bezier curves 
if we can show that for every t G [0,1] the interpolation (^ 1 ,^ 2 ) ^ di[u\^U2^f) 
exists and is a well-defined, continuous map from ( L 2 (D )) 2 to L 2 (D). To keep the 
exposition compact, we argue here formally and do not expand the arguments for the 
proper weak notion of the material derivative. 

Trouve and Younes [TY05al Theorem 6 ] have shown the existence of a geodesic 
between any U\,U2 G L 2 (D ), so for the well-definedness of X(ui,U2,t) we only need 
to show the well-definedness of evaluating such a geodesic at any time t G [0,1]. From 
HTY05all we know that any geodesic zz(t) te [ 0 ,i] is associated with a family of 

diffeomorphisms and the generating motion field v^(t) = and u(t) te [ 0 ,i] 

and 0 , 1 ] minimize the energy 


£[VV)te[o,i]> u Wte[o,i]] = 



L[v^(t),v^(t)] + da: dt, 
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i.e., £[u(t) te[ 0 j i]] = mm^ Wte[0 £[ip(t) te[0A] ,u(t) te[0tl] }. As in IBER15I , forgiven 
u(0),'«(!), and ^(f) te [ 0 j i], from BTY05al Theorem4] and ITY05al Theorem2] one ob- 
tains an explicit representation of the minimizer u(t) = o),u(i) (t) of £ [ / 0 (t) te [ 0 , 1 ] > •] 

with 


^"0,u(O),n(l) (f) 


u(0) + (ix(l) O ^(1) - tx(0)) 


Jo (det Dip) 1 (s)ds 
/ 0 ' (det Dip)~ 1 {s) ds 


o#) ^ 


Thus, taking i/(0) = xxi, xx(l) = ^ 2 , and ^ the diffeomorphism family associated with 
the geodesic, the evaluation of the geodesic at t is given by u^ jUljU2 (t) G L 2 (D). 

As for the continuity of (ui, U 2 ) T(u\,U 2 , t), consider a sequence (ui jk , ^ 2 ,/c) 

converging to ( 1 x 1 , 1 / 2 ) in ( L 2 (D )) 2 . Let us express the path energy in terms of the 
diffeomorphisms according to 

V>(*)te[o,i] ^ f fc [^Wte[o,i]] : = ,k,U2,k (f)te[o,i]] 5 

^(*)te[o,i] ^ 0,1]] •= ^[^Wte[0,i],^,ni,n 2 (f)te[0,i]] • 

It is relatively straightforward to show that minimizers 0 , 1 ] °f £ k converge 

against minimizers (V ; oo(f))te[o,i] of £°°, arguing by T-convergence of to £°° and 
equicoerciveness of the £ k . Note that the paths u^ kjUl k ^ U2 k and of images 

are geodesics. As shown in |TY05a], the ^fc(t)te[o,i] an d their inverses are even uni¬ 
formly bounded in C °’^ ([0,1], C 1,a (D)) for every a < m — 1 — | and thus admit 
strongly convergent subsequences in C 0,/3 ([0,1], C' 1 ,Q! (D)) for every 0 < (3 < \ and 
a < m — 1 — | (which we again index by k for simplicity). This strong convergence 
of the ipk and their inverses now implies the desired convergence u^ kjUl fc;U2 k (t) 
u ^oo,ui,u 2 (t) = Z{ u u u 2 ,t) in L 2 (D). Indeed, this follows from the definition of the 
U ^k,u lik ,u 2i h since for an y u in L 2 (D ) and (j) k </), (\) k X -> </> _1 in we 

have \\u k o^ k -uo(j )\\ L 2 < || det D^ 1 \\ L oo \\u k -u \\ L 2 + \\uocj) k -uo(f)\\ L 2 . The sec¬ 
ond integral converges to 0 because of the convergence of </>& in C 0 ,/ 3 ((0, l),C 1,a (D)). 
□ □ 


3.2 Discrete Bezier curves 

In analogy to discrete geodesics, we define a discrete Bezier curve 

(B K (u° 0 ,...,u° n ,k)) k =o,..., K 


by replacing the continuous operations in the defining algorithm by discrete ones. In 
explicit, a discrete iG-Bezier curve of degree n is a discrete path of images (uo ,..., uk) 
defined by n + 1 control points Uq,...,u^. For the evaluation we replace the continu¬ 
ous geodesic interpolation X(u, u , t) by the discrete geodesic interpolation I K (u, u , k) 
with t = and obtain the recursive relation 

B K (ui,... ,uj,k ) = I K (B K (ui, ..., uj-i, k), B K (u i+1 ,.. .,uj,k), k). 

For the actual computation we use the hierarchical algorithm (cf. Fig. [2]). 

If we restrict ourselves to the space of images in L 2 (D ) which are continuous up 
to a null set, we also obtain in the discrete context an existence and continuity result. 
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Figure 2: The Discrete de Casteljau algorithm with a schematic sketch for n — 3, 
K = 7, and k = 4 (j = 1 green, j = 2 orange, j = 3 blue). 


Theorem 3.2 (Existence and stability of discrete Bezier curves). Let D be a Lipschitz 
domain and n, K > 1. For any n + 1 images ..., G L 2 (D) there exists 
a discrete Bezier curve ~B K (u ®,..., fc) G L 2 (D), k = 0,...,iC. Furthermore , if 
(uqj, ..., )j—1,2,... ^ ^ sequence of control points converging against (rz[], ..., ^) 

in (L 2 (D)) n+1 , exist Bezier curves B k (izqj, ..., j, •) of which a subse¬ 
quence converges pointwise against a Bezier curve ~B k ..., •). 

Proof The existence follows directly from the existence result IBER15I Theorem 3.4] 
for discrete geodesics. The proof of continuity follows the line of argumentation in the 
proof of Theorem |3T] using compactness for the underlying family of deformations 
and the fact that for two images uo , uk and a discrete connecting geodesic with associ¬ 
ated deformations </>i,..., fx, the images = l K (uo, uk , fc) along the geodesic can 
be expressed as the unique solution of a block tridiagonal system of operator equations 
1BER15I Proposition 3.2] whose entries continuously depend on fi ,..., fx- n □ 

Let us remark that the convergence of discrete Bezier curves for K —>> oo to a con¬ 
tinuous Bezier curve for the same set of input shapes is still open. In fact, in 1BER151 
T-convergence of discrete geodesics is proven in the L 2 ((0,1) x D) topology, which 
is too weak to control the convergence of the interpolated images in the de Casteljau 
algorithm pointwise with respect to time. 

To compute discrete geodesics numerically we have to introduce a suitable space 
discretization. To this end, we use a finite element discretization both for the family of 
images and for the family of deformations. For details we refer to 1BER151 . 
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Figure 3: Piecewise discrete geodesic (middle) and cubic Bezier curve (bottom). K = 
8, 5 = 5- IQ- 2 , 7 = 1(T 3 . 
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Figure 4: Piecewise discrete geodesic (middle) and cubic Bezier curve (bottom). K = 

8,6 = 5 - 1CT 3 , 7 = icr 3 . 

Fig. [5] and Fig. [4] show discrete Bezier curves for two different test cases to 
highlight some of the general characteristics of Bezier curves in the space of im¬ 
ages equipped with the Riemannian structure of the metamorphosis approach. In both 
cases the initial image at t = 0 and the final image at t = 1 coincide with the cor¬ 
responding control images, where as for the other control images the cubic Bezier 
curves approximately recovers them at times t = | and t = |. The Bezier curve 
is smooth in time and structures in the images are smoothly blended and transported 
(cf. also the accompanying video sequence). The impact of the control images is 
global in time. Thus particular features of a single control image persist over the 
whole time interval (0,1). In Fig. [3] this can be best observed for the rectangu¬ 
lar outline of the second and third control images and for the lower tip of the tri¬ 
angle. In Fig. [4] one obtains a simultaneous rotation of the emerging grey ellip¬ 
soid and a continuous fading in and out of the white and black circle, respectively. 
Compared to this the piecewise geodesic 
interpolation is not smooth in time and 
local features of the input images spread 
just over two consecutive time inter¬ 
vals. Note that the metamorphosis 
Bezier curve in Figure[3] is not solely 
generated by a transport. Indeed, the generation of the lower tip of the triangle is 
not possible via a flow of diffeomorphisms. Hence, the triangle tip is generated via 
intensity modulation (cf. Fig. [5]). 



Figure 5: Zoom into marked regions in Fig¬ 
ure [3 ] 


4 Applications 

In this section we present two applications of Bezier curves in the space of images. 
First, we use our approach for the modulation of human face interpolation. Then we 
compute discrete Bezier curves to reconstruct animations of shapes. 
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Figure 6 : Piecewise discrete geodesic (middle) and quadratic Bezier curve (bottom) 
between human faces. K = 8 , S = 7.5 • 10 -3 , 7 = 10 -3 . 

4.1 Modulation of Interpolation Path 

Bezier curves can be used to modulate a (geodesic) interpolation path between two 
given images via prescribing additional control images, which are used to let additional 
image pattern appear and spread along a smooth and natural looking path connecting 
the two end images. Figures [ 6 ] and [T] display a quadratic and a cubic Bezier curve be¬ 
tween different portraits. Furthermore the Bezier curves are compared with piecewise 
geodesic interpolations between the control images. The resulting path between the 
start and the end image ((A,C) and (A,D), respectively) consists of natural looking fa¬ 
cial image for all times. The resulting path in image space is smooth in time. Features 
of the control images are clearly visible and significantly spread in time over the image 
sequence. E.g .in Fig. [ 6 ] the hairstyle of face C is clearly visible in the images at time 
steps 2 — 7 of the Bezier curve. Furthermore, the shape of the chin of face A carries 
over up to the time step 5 — 6 and the moustache from face B spreads over the whole 
sequence from time step 1 to 7. 

4.2 Image Sketches as Control Shapes for Animation 

Discrete Bezier curves can also be used to design animation paths in the space of im¬ 
ages based on very moderate user interaction. In fact, to animate an object in an image 
the user might cut the object into pieces and generate 
control images via translation, rotation or scaling of 
the different pieces. Then, de Casteljau’s algorithm 
is used to generate from these sketches a smooth an¬ 
imation path in the space of images, which approx¬ 
imately recovers the poses described by the control 
images. As a proof of concept we consider in Fig. 

[ 8 ] four control images (A-D), which are generated 
from a very simple model of an arm (top left). While 
a geodesic curve between the first (A) and the last 
image (D) prefers to collapse the forearm (second 
row), the additional control points (B,C) give further 
guidance to reconstruct the bending of the arm along 



Figure 7: Modeling of control 
images. 
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Figure 8: A simple model of an arm (top left) is split up into three pieces, they are 
reconfigured into four control images (A-D). A geodesic interpolation of the images 
(A) and (D) (middle row) and a discrete cubic Bezier curve with K = 8, S = 10 -2 , 
7 = 10 -3 is computed (bottom row). 
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Figure 9: An animation sequence based on the image of an elephant (www.fotolia.com 
©eyetronic) via discrete cubic Bezier curve is shown for K = 8, S = 4 • 10 -3 , 
7 = 10 -3 . The underlying control images are depicted in Fig. [7] 


a discrete cubic Bezier curve (third row). Next, we consider in Fig. [7] a photograph of 
an elephant cut into three pieces (A), which are reconfigured into four control images 
(A-D). These control images are then taken into account to compute a discrete cubic 
Bezier curve (Figure [9]). The resulting sequence represents the turning and stretching 
of the head and the raising of the forelegs. 


5 Conclusions 

We have defined Bezier curves on the image manifold equipped with the metamorpho¬ 
sis metric. Based on the notion of a discrete metamorphosis path energy and a cor¬ 
responding discrete geodesic interpolation, de Casteljau’s algorithm allows the robust 
and efficient computation of discrete Bezier curves in the space of images. Striking 
features of this approach are the smoothness in time and the global impact of features 
of the control images. This approach can be regarded as a conceptual study for curve 
modeling on shape spaces. Indeed, Bezier curves are just one prominent example for 
classes of curves constructed via hierarchical convex combination along straight lines. 
Perspective future generalizations include the Neville-Aitken scheme for general poly¬ 
nomials and B-Spline curves. 
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